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I.  INTRODUCTION 

The  Incomplete  Beta  Function,  Ix(a,  b)  [1,  p.  263,  p.  944],  is  defined  by 

Ix(a,  b)  =  G(a,  b)  /  ta— 1  (1  —  t)b  1  dt,  a  >  0,  b  >  0,  0  <  x  <  1  (1) 

Jo 

B(a,  b)  =  1/G(a,  b)  =  [  ta_1  (1  -  t)15"1  dt  =  T(a)  r(b)/r(a  +  b),  (2) 

Jo 

where  the  gamma  function  r(u)  is  given  by 

/»oo 

r(u)  =  /  e_ttu_1dt,  u  >  0,  [1,  p.  255].  (3) 

Jo 

The  complement  of  Ix(a,  b)  is  given  by 

Iy(b,a)  =  1  — Ix(a,b),  y  =  1  -  x.  (4) 

Ix(a,  b)  is  numerically  evaluated  by  the  Fortran  95  subroutine  called  DBRAT  contained 
in  the  NSWC  Mathematics  Library  [7,  p.88].  It  requires  a, b,x, y  as  input,  where  both  x 
and  y  are  specified,  so  that  a  value  for  either  x  or  y  that  is  smaller  than  e  can  be  used. 
Here  e  is  the  smallest  positive  double-precision  number  for  the  computer  in  use  such  that 
1  +  e  >  0  (see  [2],  [3]).  In  our  use  of  an  IBM  PC,  e  =  2~52  ~  2.220446E(— 16). 

Ix(a,  b)  appears  in  many  branches  of  science,  including  atomic  physics,  fluid  dynamics, 
transmission  theory,  lattice  theory,  and  operations  research  [4],  It  is  perhaps  best  known 
in  statistics  by  its  direct  connection  to  the  F  —  Distribution,  P(F0|^i,  z/2)  [1,  p.946].  The 
F  —  Distribution  and  its  complement  are  given  by 

P(Foh,  v2)  =  zql/2  vv22'2  G(iq/2,  iV2)  f  °  F(i/1“2)/2(z/2  +  z/i  F)~^+^2  dF,  (5) 

Jo 

Q(F0|z/i,  v2)  =  1  -  P(F0|z/i,  u2).  (6) 

Let 

=  2a,  v2  =  2b,  F  =  b(l  —  t)/(at),  F0  =  (by)/(ax).  (7) 

Then,  it  easily  follows  that  x  =  b/(b  +  aF0),  y  =  (aFo)/(b  +  aF0)  and 

P(F0|u,  ^2)  =  Iy(a,b),  Q(F0|i*i,  v 2)  =  Ix(b,  a).  (8) 

With  a  =  1/2,  v  —  v2,  the  F— Distribution  reduces  to  the  Student's  t  —  Distribution, 

( A(t0 1 zv)  [1, p.948].  Thus 

P(|t|  <  to)  =  A(toW  =  5)  J\l  +  t2nH"+1),2dt  =  Iy(l/2,b),  (9) 
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where 

to=V^,  y  ~  tg/ [y  +  tg).  (10) 

In  addition,  Ix(a,  b)  is  also  related  to  the  Binomial  Distribution,  E  [1,  p.  960],  namely 

E(n,r,x)  =  ^2  _x)n~1  =  Ix(r,n-r  +  1).  (11) 

i=r  '  ' 

The  objective  in  this  report  is  to  give  a  numerical  algorithm  to  obtain  the  inverse  of  Ix, 
that  is,  to  find  x  (and  y)  given  a,  b,  Ix,  Iy.  In  addition,  a  Fortran  95  computer  program  will 
be  described  that  determines  the  smaller  of  x  and  y  within  a  prescribed  number  of  significant 
digits.  For  reasons  that  required  the  input  of  both  x  and  y  for  DBRAT  [2],  we  will,  similarly, 
now  require  both  Ix  and  Iy  as  input. 

II.  ALGORITHM  FOR  x  AND  y 

In  order  to  obtain  values  of  x  and  y  with  a  prespecified  relative  accuracy,  the  smaller  of 
the  two  is  always  the  one  computed,  and  the  other  is  determined  as  its  one’s  complement. 
Also,  in  referring  to  Ix(a,  b)  or  Iy(b,  a),  it  is  usually  the  smaller  that  is  used.  Note  that 
Ix(a,  b)  <  Iy (b ,  a)  does  not  necessarily  imply  that  x  <  y.  See  the  example  directly  following 
(18).  Therefore  let  Izn  (Izo)  denote  the  new  (previous)  iterate  for  Iz,  where  z  =  x  or  y  or  xy, 
with 


xy  =  mini 
Also  let 

A,y),  Ixy  =  min(Ix,  Iy), 

Ixyn 

=  min(Ixn,  Iyn)  Ixyo  = 

1— 1 

O 

X 

1— 1 

£ 

a 

(12) 

Previous  iterate  for  x  = 

xo. 

New  iterate  for  x  = 

xn. 

Previous  iterate  for  y  = 

yo. 

New  iterate  for  y  = 

yn. 

(13) 

Previous  iterate  for  xy  = 

xyo. 

New  iterate  for  xy  = 

xyn. 

The  numerical  evaluation  of  xy  from  Ixy  for  given  a  and  b  is  essentially  an  iterative 
rootfinding  process.  For  this  purpose,  the  Newton-Raphson  (N-R)  procedure  [5,  p.  129] 
worked  very  well.  By  this  we  mean  that  it  was  successfully  used  to  find  xy  numerically  over 
a  large  range  of  input.  However,  before  one  can  use  N-R  with  assurance,  a  so-called  “domain 
of  attraction”  (DOA)  must  be  determined,  which  contains  an  estimate  for  xy  from  which  a 
few  N-R  iterations  will  converge  to  an  acceptable  approximation  for  xy.  Extensive  testing 
determined  that  if  Ixyn  satisfies 

|Ixyn  Ixy  |  epsl  *  Ixy;  epsl  1(  2)  ,  (14) 

then  xyn  is  in  DOA.  With  three  independent  variables,  a,  b,  Ixy,  this  is  the  difficult  phase  of 
the  analysis. 

1A(— B)  =  A  *  10— B,  A  and  B  positive  numbers. 


2 


NSWCDD/TR— 05/91 


For  our  purposes  it  is  useful  to  always  have  an  upper,  xyH,  and  lower,  xyL,  bound  for  xy 
such  that 

xyL  <  xy  <  xyH.  (15) 

Initially,  we  start  with  xyL  =  0  and  xyH  =  1.  The  first  approximation  for  xy  is  given  by 
either 


xo  =  [aIxB(a,  b)]ly/a, 

a  <  b, 

(16) 

yo  =  [b  Iy  B(a,  b)]1/b, 

b  <  a. 

(17) 

Estimate  (16)  is  particularly  good  for  small  a,  and  (IT)  for  small  b.  The  heuristic  argument 
for  (16)  is  based  on  the  relationship  [1,  p.944] 

Ix(a,b)  =  G(a,b)xa(l  -  x)b  +  Ix(a  +  l,b).  (18) 

cl 

Since  Ix(a,  b)  — >  1  as  a  — >  0,  then  for  small  a  and  reasonable  values  of  Ix,  x  will  be  small  so 
that  the  last  term  in  (18)  is  negligible  relative  to  the  first,  thus  yielding  (16).  For  example, 
given  Ix  =  .9  and  a  =  .01,  b  =  2  (from(16)),  xo  ~  9.80047(— 6)  with  x  =  9.801426(— 6).  An 
argument  similar  to  the  one  above  leads  to  (17). 

Improved  estimates  for  xy  are  attempted  by  the  following:  If 

Ixyo  <  Ixy  (see  (13)),  (19) 

then  xyn  =  2*xyo  <  1.  This  doubling  procedure  is  continued  for  a  maximum  of  12  cycles  as 
long  as  (19)  is  satisfied  with  a  resulting  improvement  in  xyL.  If  (19)  is  not  satisfied  initially, 
then  xyn  =  xyo/2.  This  halving  procedure  is  continued  for  a  maximum  of  12  cycles  as  long 
as  (19)  is  not  satisfied  with  a  resulting  improvement  in  xyH. 

At  this  stage  xyn  is  assigned, 

xyn  =  (xyH  +  xyL)/2,  (20) 

so  that  an  improved  value  of  xyH  or  xyL  is  found  depending  on  whether  Ixyn  is  greater  or 
less  than  Ixy.  This  procedure  is  cycled  using  (20)  a  maximum  of  35  times  unless,  at  some 
stage,  (14)  is  satisfied  or 

|xyH  —  xyL|  <  eps2  *  xyn,  eps2  =  1(— 13).  (21) 

If  (21)  holds,  then  xyn  is  accepted  as  a  satisfactory  result  for  xy. 

If  only  (14)  holds,  then  xyn  is  in  DOA  and  the  N-R  procedure  is  called.  Setting  xyn  to 
xyo  (see  (13)),  the  new  estimate  for  xy  is  computed  from  (N-R)  using  the  following  algorithm 
(see  next  page). 
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Let 


Now: 


D 


xyo 


xoa  (1  —  xo)b/ ( B(a,  b)  xo  yo )  for  xo  <  yo 
yob  (1  —  yo )a/ ( B(a,  b)  xo  yo )  for  yo  <  xo. 


If  Ix  <  Iy  then 

w  =  (Ixyo  -  Ix)/Dxyo, 
If  xo  <  yo  then 

xyn  =  xn  =  xo  —  w 
Else 

xyn  =  yn  =  yo  +  w 
Enclif 

Else  (Iy  <  Ix) 

W  =  (Ixyo  Iy)/DXyo) 
If  yo  <  xo  then 

xyn  =  yn  =  yo  —  w 
Else 

xyn  =  xn  =  xo  +  w 

Enclif 

Enclif 


The  N-R  procedure  is  cycled  until 


(22) 


(23) 


|  w  |  <  eps  *  xyn,  (24) 

where  eps  is  assigned  by  the  user.  An  eps  =  5(— 10)  requires,  for  convergence,  no  more  than 
4  N-R  iterations  for  a  and/or  b  as  large  as  1(8)  and  Ixy  as  small  as  1(-10). 

III.  COMPUTER  PROGRAM  FOR  x  and  y— INVBETA  (  a,  b,  Ix,  Iy,  x,  y,  eps,  Ind) 

INVBETA  (a,  b,  Ix,  Iy,  x,  y,  eps,  Ind )  denotes  the  Fortran  95  [6]  subroutine  with  its  call 
line.  All  variables  in  the  call  line  are  specified  double-precision,  except  the  integer  Inch  The 
inputs  are  a,  b,  Ix,  Iy  and  eps.  These  have  been  defined  in  the  previous  sections.  The  outputs 
are  x,  y  and  Inch  The  input  eps  assigns  the  accuracy  desired  in  x  and  y.  If  a  relative  error  less 
than  1(-N)  is  desired  in  x  and  y,  then  eps  =  5(-N-l).  The  value  of  N  is  limited  to  a  positive 
integer  no  larger  than  10.  The  integer  output  Ind  specifies  whether  the  output  x  and  y  are 
acceptable.  The  tablulation  below  indicates  the  values  Ind  can  take  and  their  meanings. 

Ind  =  0,  Results  for  x  and  y  are  acceptable. 

Ind  =  1,  a  <  0  and/or  b  <  0  are  unacceptable,  a  >  0,  b  >  0  required. 

Ind  =  2,  Ix  +  Iy  ^  1  is  unacceptable,  Ix  +  Iy  =  1  required. 

Ind  =  3,  Ix  <  0  or  Ix  >  1  are  not  acceptable. 


4 


NSWCDD/TR— 05/91 


Ind  =  4,  Iy  <  0  or  Iy  >  1  are  not  acceptable. 

Ind  =  100,  x  is  estimated  to  be  less  that  l(-200). 

Ind  =  —100,  y  is  estimated  to  be  less  that  l(-200). 

Ind  =  10,  More  than  20  iterations  of  N-R  required  on  x,  returns  to  halving  stage. 

Ind  =  11,  More  than  20  iterations  of  N-R  required  on  y,  returns  to  halving  stage. 
The  Fortran  source  hie  IBETA.FOR  contains  INVBETA  and  43  supporting  routines  from 

[7]- 

IV.  NUMERICAL  RESULTS 

In  Table  1  below,  the  Erst  four  columns  contain  inputs  for  INVBETA.  The  next  two 
columns  contain  outputs  x,  y,  and  the  last  column  contains  the  total  number  of  calls  to 
DBRAT.  A  row  of  data  refers  to  one  case,  or  call  to  INVBETA.  There  are  23  cases  listed. 
The  accuracy  parameter  was  set  at  eps  =  5(— 11). 


TABLE  1.  NUMERICAL  RESULTS 


l(-2) 

10 

.99 

l(-2) 

,2748462647(-l) 

.9725153735 

7 

10 

.10 

.9999 

l(-4) 

1.0 

.6356342847(-41) 

3 

2 

5 

.1 

.9 

.9259525891(-1) 

.9074047411 

11 

2 

5 

.5 

.5 

.2644499833 

.7355500167 

8 

2 

5 

.9 

.1 

.5103163066 

.4896836934 

9 

2 

5 

.999 

l(-3) 

.8186138669 

.1813861331 

9 

10 

25 

.1 

.9 

.1914227060 

.8085772940 

8 

10 

25 

.5 

.5 

.2815933420 

.7184066580 

10 

10 

25 

.9 

.1 

.3854635784 

.6145364216 

12 

10 

25 

.999 

l(-3) 

.5413352819 

.4586647181 

14 

50 

125 

.1 

.9 

.2425526829 

.7574473171 

13 

50 

125 

.5 

.5 

.2848963507 

.7151036493 

12 

50 

125 

.9 

.1 

.3299331965 

.6700668035 

15 

50 

125 

.999 

l(-3) 

.3969077527 

.6030922473 

11 

500 

925 

.1 

.9 

.3347225032 

.6652774968 

14 

500 

925 

.5 

.5 

.3508074095 

.6491925905 

12 

500 

925 

.9 

.1 

.3671216035 

.6328783965 

15 

500 

925 

.999 

l(-3) 

.3904834123 

.6095165877 

16 

50000 

2500000 

.1 

.9 

,1949665343(-1) 

.9805033466 

19 

50000 

2500000 

.5 

.5 

.1960771754(-1) 

.9803922825 

14 

50000 

2500000 

.9 

.1 

.1971919420(-1) 

.9802808058 

17 

50000 

2500000 

.999 

l(-3) 

.1987722641(-1) 

.9801227736 

20 

50000 

2500000 

.9999999999 

i(-io) 

.2016513068(-1) 

0.9798348693 

21 
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OPERATIONS  RESEARCH 
VIRGINIA  COMMONWEALTH  UNIVERSITY 
PO  BOX  842014 
RICHMOND  VA  23284-2014 

DR  JEFFREY  BIRCH  1 

DEPT  OF  STATISTICS 

VIRGINIA  TECH 

4 17- A  HUTCHESON  HALL 

BLACKSBURG  VA  24061 

DR  DONALD  RICHARDS  1 

DEPT  OF  STATISTICS 
UNIVERSITY  OF  VIRGINIA 

103  HALSEY  HALL  1 

CHARLOTTESVILLE  VA  22904-4135 

DR  ED  CARLSTEIN  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  NORTH  CAROLINA— 

CHAPEL  HILL 

NEW  WEST  BUILDING 

CHAPEL  HILL  NC  27599-3260 

DR  MARCIA  GUMPERTZ  1 

DEPT  OF  STATISTICS 
NORTH  CAROLINA  STATE  UNIVERSITY 
RALEIGH  NC  27695-8203 

DR  JAMES  D  LYNCH  1 

DEPT  OF  STATISTICS 
UNIVERSITY  OF  SOUTH  CAROLINA 
COLUMBIA  SC  29208 

DR  JAMES  R  SCHOTT  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  CENTRAL  FLORIDA 

PO  BOX  162370 

ORLANDO  FL  32816-2370 

DR  MYLES  HOLLANDER  1 

DEPT  OF  STATISTICS 

FLORIDA  STATE  UNIVERSITY 

PO  BOX  118545 

GAINESVILLE  FL  32611-8545 


DR  JACK  D  TUBBS  1 

INSTITUTE  OF  STATISTICS 
BAYLOR  UNIVERSITY 
PO  BOX  97140 
WACO  TX  76798 
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DR  JAMES  STAPLETON  1 

DEPT  OF  STATISTICS  AND  PROBABILITY 
MICHIGAN  STATE  UNIVERSITY 
EAST  LANSING  MI  48824 

DR  ENSOR  KATHERINE  1 

DEPT  OF  STATISTICS 

RICE  UNIVERSITY 

P.O.  BOX  1892  MS  138 

HOUSTON  TX  77251-1892 

MS  BETSEY  COGSWELL  1 

DEPT  OF  STATISTICS 
HARVARD  UNIVERSITY 
1  OXFORD  ST 
CAMBRIDGE  MA  02138-2901 

DR  EDWARD  R  SCHEINERMAN  1 

DEPT  OF  MATHEMATICS 

JOHNS  HOPKINS  UNIVERSITY 

104  WHITEHEAD  HALL 

BALTIMORE  MD  21218 

DR  BENJAMIN  KADEM  1 

DEPT  OF  MATHEMATICS 
UNIVERSITY  OF  MARYLAND 
COLLEGE  PARK  MD  20742 

DR  WAYNE  A  MUTH  1 

DEPT  OF  STATISTICS  AND  COMPUTER  SCIENCE 
WEST  VIRGINIA  UNIVERSITY 
PO  BOX  6330 

MORGANTOWN  WV  26506-6330 

MS  SHARON  DINGESS  1 

DEPT  OF  STATISTICS 

CASE  WESTERN  RESERVE  UNIVERSITY 

10900  EUCLID  AVE  YOST  323 

CLEVELAND  OH  44106-7054 

DR  DOUGLAS  WOLFE  1 

DEPT  OF  STATISTICS 

OHIO  STATE  UNIVERSITY 

1958  NEIL  AVE 

COLUMBUS  OH  43210-1247 

DR  ERIK  V  NORDHEIM  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  WISCONSIN-MADISON 

1210  W  DAYTON  ST 

MADISON  WI  53706-1613 

DR  GARY  OEHLERT  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  MINNESOTA 

224  CHURCH  ST  SE  STE  313 

MINNEAPOLIS  MN  55455-0493 


MS  EDITH  E  LANDIN  1 

DEPT  OF  STATISTICS 

IOWA  STATE  UNIVERSITY 

102  SNEDECOR  HALL 

AMES  IA  50011-1210 

DR  DALE  ZIMMERMAN  1 

DEPT  OF  STATISTICS  AND  ACTUARIAL  SCIENCE 
UNIVERSITY  OF  IOWA 
241  SCHAEFFER  HALL 
IOWA  CITY  IA  52242-1409 

DR  JOHN  E  BOYER  1 

DEPT  OF  STATISTICS 

KANSAS  STATE  UNIVERSITY 

DICKENS  HALL  101 

MANHATTAN  KS  66506 

DR  DHARAM  V  CHOPRA  1 

DEPT  OF  STATISTICS 

WICHITA  STATE  UNIVERSITY 

1845  N  FAIRMOUNT 

WICHITA  KS  67260-0033 

MS  CAROLYN  J  COOK  1 

DEPT  OF  STATISTICS 

COLORADO  STATE  UNIVERSITY 

101  STATISTICS  BLDG 

FORT  COLLINS  CO  80523-1877 

DR  DENNIS  L  YOUNG  1 

DEPT  OF  MATHEMATICS  AND  STATISTICS 
ARIZONA  STATE  UNIVERSITY 
TEMPE  AZ  85287-1804 

DR  DOOD  KALICHARAN  1 

DEPT  OF  STATISTICS 

COLUMBIA  UNIVERSITY 

2990  BROADWAY  MC  4403 

NEW  YORK  NY  10027 

DR  LU  ANN  CUSTER  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  MICHIGAN 

439  WEST  HALL 

ANN  ARBOR  MI  48109-1092 

DR  ROBERT  SMYTHE  1 

DEPT  OF  STATISTICS 

OREGON  STATE  UNIVERSITY 

KIDDER  HALL  44 

CORVALLIS  OR  97331 
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DR  KENT  ESKRIDGE  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  NEBRASKA- LINCOLN 

926  OLDFATHER  HALL 

LINCOLN  NE  68588 

DR  WILLIAM  D  WARDE  1 

DEPT  OF  STATISTICS 

OKLAHOMA  STATE  UNIVERSITY 

301  MSCS  BLDG 

STILLWATER  OK  74078-1056 

DR  ISHWAR  V  BASAWA  1 

DEPT  OF  STATISTICS 
UNIVERSITY  OF  GEORGIA 
ATHENS  GA  30602 

DR  WILLIAM  GRIFFITH  1 

DEPT  OF  STATISTICS 
UNIVERSITY  OF  KENTUCKY 
LEXINGTON  KY  40506-0027 

DR  REBECCA  W  DOERGE  1 

DEPT  OF  STATISTICS 

PURDUE  UNIVERSITY 

1399  MATHEMATICAL  SCIENCE  BLDG 

WEST  LAFAYETTE  IN  47907 

DR  LARRY  WASSERMAN  1 

DEPT  OF  STATISTICS 
CARNEGIE  MELLON  UNIVERSITY 
PITTSBURGH  PA  15213 

DR  SATISH  IYENGAR  1 

DEPT  OF  STATISTICS 
UNIVERSITY  OF  PITTSBURGH 
2703  CATHEDRAL  OF  LEARNING 
PITTSBURGH  PA  15260 

DR  DIRK  MOORE  1 

DEPT  OF  STATISTICS 

TEMPLE  UNIVERSITY 

SPEAKMAN  HALL  006-00 

PHILADELPHIA  PA  19122 

DR  DIPAK  K  DEY  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  CONNECTICUT 

196  AUDITORIUM  ROAD 

STORRS  CT  06269 


DR  ANDREW  BARRON  1 

DEPT  OF  STATISTICS 

YALE  UNIVERSITY 

PO  BOX  208290 

NEW  HAVEN  CT  06520-8290 

DR  KESAR  SINGH  1 

DEPT  OF  STATISTICS 

RUTGERS  UNIVERSITY  OF  NEW  JERSEY 
501  HILL  CENTER  BUSCH  CAMPUS 
PISCATAWAY  NJ  08854 

DR  ROBERT  SMYTHE  1 

DEPT  OF  STATISTICS 
GEORGE  WASHINGTON  UNIVERSITY 
WASHINGTON  DC  20052 

DR  STEPHEN  L  BIEBER  1 

DEPT  OF  STATISTICS 

UNIVERSITY  OF  WYOMING 

PO  BOX  3332 

LARAMIE  WY  82071-3332 

DR  HOWARD  B  CHRISTENSEN  1 

DEPT  OF  STATISTICS 
BRIGHAM  YOUNG  UNIVERSITY 


230A  TMCB 
PROVO  UT  84602 

DR  ELGIN  PERRY  1 

2000  KING’S  LANDING  RD 
HUNTINGTOWN  MD  20639-9743 

INTERNAL 

B60  TECHNICAL  LIBRARY  3 

B10  DR  ALAN  BERGER  1 

B10  DR  JOHN  CRIGLER  1 

B10  DR  WILLIAM  FARR  1 

B10  NGA  PHAM  1 

G24  DR  THOMAS  GOSWICK  1 

G33  CHARLES  GARNETT  1 

G33  ORGAL  HOLLAND  1 

T30  ROBERT  HILL  1 

T30  DR  ARMIDO  DIDONATO  2 

T31  DR  JEFFREY  BLANTON  1 

T32  DR  ROBERT  MCDEVITT  1 

T41  DR  MICHAEL  RUDZINSKY  1 

T51  WILLIAM  ORMSBY  1 

T54  DAVID  CLAWSON  1 

T505  STEVEN  ANDERSON  1 
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